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Abstract 

The Conjugate Gradient method (CGM) is known to be the fastest generic iterative 
method for solving linear systems with symmetric sign definite matrices. In this paper, we 
modify this method so that it could find fundamental solitary waves of nonlinear Hamil- 
tonian equations. The main obstacle that such a modified CGM overcomes is that the 
operator of the equation linearized about a solitary wave is not sign definite. Instead, it has 
a finite number of eigenvalues on the opposite side of zero than the rest of its spectrum. We 
present versions of the modified CGM that can find solitary waves with prescribed values of 
either the propagation constant or power. We also extend these methods to handle multi- 
component nonlinear wave equations. Convergence conditions of the proposed methods are 
given, and their practical implications are discussed. We demonstrate that our modified 
CGMs converge much faster than, say Petviashvili's or similar methods, especially when 
the latter converge slowly. 
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1 Introduction 



Solitary waves of most nonlinear evolution equations can be found only numerically. In one 
spatial dimension, a variety of numerical methods for finding solitary waves are available, 
with the shooting and Newton's methods being probably the most widely used ones. How- 
ever, in two and three spatial dimensions, the shooting method is no longer applicable, and 
the complexity of programming of Newton's method increases considerably compared to the 
one-dimensional case. Also, the accuracy of Newton's method is only algebraic in the size 
of the spatial discretization step Ax (typically, it is O ((Ax) 2 )). This accuracy may be too 
low for some applications, since in two or three spatial dimensions, one tends to take larger 
discretization steps than in one dimension in order to limit the computational time and the 
size of stored numeric arrays. Therefore, it is desirable to have a method that would: (i) 
be easier to program than Newton's method; (ii) have the exponential accuracy in Ax, as 
spectral methods; and (iii) converge sufficiently fast under known conditions. 

A number of such methods are available. The well-known iterative method proposed by 
Petviashvili pQ can be used to find stationary solitary waves of equations with power-law 
nonlinear ity: 

-Mu + u p = 0, u(x) -> as |x| oo, (1.1) 

where u is the real- valued field of the solitary wave, M is a positive definite and self-adjoint 
differential operator with constant coefficients, and p is a constant. For example, the solitary 
wave of the nonlinear Schrodinger equation in D spatial dimensions, 

iU t + V 2 U +\U\ 2 U = 0, U(\x\ -^oo) ->0, 

y2 = ^ (1-2) 
dx\ dx 2 D ' 

upon the substitution U(x, t) = e* M *u(x) reduces to Eq. (11.11) with p = 3 and 

M = /i-V 2 . (1.3) 

The parameter [i > is referred to as the propagation constant of the solitary wave. 
Convergence conditions of the Petviashvili method for Eq. (jl.ip were established in [2]. 

Studies of solitary waves in nonlinear photonic lattices and Bose-Einstein condensates 
motivated the interest in finding solitary waves of equations that have a more general form 
than (ITT]) : 

- Mu + F(u,x) = 0, it(x) -> as |x| -> oo, (1.4) 

where F(u, x) is a real function. For example, solitary waves in nonlinear photonic lattices 
with Kerr nonlinearity satisfy the following equation of the form (11,51k 

V 2 u + Vb(cos 2 x + cos 2 y) u + u 3 = \iu , (1-5) 

where the second term accounts for the effect of the periodic potential provided by the 
lattice. Various modifications of the Petviashvili method for obtaining solutions of (jl.4p 
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were proposed; see, e.g., [3], 01 El E] . However, convergence conditions of those modifications 
were not studied. In Ref. [3 [8], J. Yang and the present author proposed two alternative 
iterative methods which satisfied conditions (i)-(iii) stated before Eq. (jl.ip ; in particular, 
the conditions of their convergence to fundamental (see below) [7] and both fundamental 
and non-fundamental solitary waves (also known as ground and excited states, respectively) 
[8] were given. In Refs. [9j[8], a technique referred there as mode elimination was proposed, 
which was shown to provide considerable acceleration of iterative methods when the later 
converged slowly to the respective solitary wave. 

In this paper, we propose yet another family of numerical methods for finding fun- 
damental solitary waves of Hamiltonian nonlinear evolution equations. Consider a linear 
system 



for the unknown vector y, where A is a real symmetric matrix and b is a known vector. Then 
the methods proposed here are extensions of the well-known Conjugate Gradient method 
(CGM) for the linear system (jl.6p to the nonlinear boundary- value problem (|l,4p . and they 
converge even faster than, e.g., the generalized Petviashvili method [7] accelerated by mode 
elimination [9]. 

The main part of this paper is organized as follows. In Section 2, we explain how 
the methods introduced in [TJ El E] for the nonlinear problem (|1.4|) are related to some 
well-known iterative methods of solving the linear system (jl.6p with a real and symmetric 
matrix A. This comparative review will lead us to explaining the main issue that needs 
to be resolved in order to extend the CGM in such a way that it could find fundamental 
solitary waves. In Section 3, we present such an extension of the CGM for the case where the 
solitary wave of the single-component Eq. (11.41) has a prescribed value of the propagation 
constant /i. In Section 4, we present a modification of the CGM for finding fundamental 
waves of (|1.4p with a prescribed value not of the propagation constant but of power, defined 
as 



where the integration is over the entire spatial domain. Let us note that another method 
for finding solitary waves with a specified value of power was previously considered in, e.g., 
[El El [121 CE] , with its convergence conditions for (II. 4p being established in [13]. We will 
refer to this method as the imaginary-time evolution method. In Section 5, we generalize 
the modified CGMs of Sections 3 and 4 to multi-component equations. In Section 6, we 
compare the performances of the generalized Petviashvili and imaginary-time evolution 
methods accelerated by the mode elimination technique [91 [8] with the performances of 
the respective modified CGMs from Sections 3-5. In Section 7 we summarize this work. 
Appendix 1 contains proofs and discussions of some auxiliary results of Sections 3 and 4. 
Appendix 2 contains an alternative method to the modified CGM of Section 3. Appendix 



Ay - b = 



(1.6) 




(1.7) 
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3 contains a sample code of a modified CGM for a two-component equation considered in 
Section 6. 

Before concluding this Introduction, we need to discuss what we will refer to as a fun- 
damental solitary wave. In fact, we are not aware of a universal definition of this term. 
In many situations, one can intuitively identify a solitary wave as being fundamental if it 
satisfies two easily-verifiable conditions. First, it is to have one main peak, with smaller 
peaks possibly existing around the main one. Second, for envelope solitary waves (e.g., those 
satisfying the complex-valued Eq. (|1.2p ). the propagation constant of the fundamental wave 
must lie in the semi-infinite spectral bandgap; see Fig.UJin Section 6. (Typically, only equa- 
tions with periodic potentials or systems describing more than two linearly coupled waves 
propagating with different group velocities have more than one spectral bandgap.) Alter- 
natively, for carrier (e.g., Korteweg-de Vries-type) solitary waves, the velocity, rather than 
the propagation constant, of a fundamental wave must lie in the semi-infinite gap. More 
rigorous characterization of fundamental waves is possible in special cases. For instance, 
when the spatial operator in the wave equation is the Laplacian (as, e.g., in (jl.2p or (jl.5p ). 
the fundamental wave is known to be nonzero for all x (see, e.g., |14|): however, for a more 
general operator as, e.g., in the Kadomtsev-Petviashvili equation, the fundamental wave 
may have zeros [15]. For some equations (see, e.g., [13 [T7]), it was rigorously shown that 
the fundamental solution minimizes a functional known as the action. 

The property of an 5-component fundamental solitary wave that we will rely on in this 
paper is that for many such waves, the linearized operator of the stationary equation has 
no more than S positive eigenvalues. To the author's knowledge, there is no general proof 
of this property even for a single-component Eq. (|1.4p . let alone for the S-component case. 
(See, however, a recent review of solitary wave stability in lattices [18].) Therefore, we sim- 
ply declare this property to be our working definition of a fundamental solitary wave for the 
purpose of this paper. Again, we do not imply that any solitary wave that has the two intu- 
itive properties described at the beginning of the previous paragraph also satisfies the above 
property about the positive eigenvalues of the corresponding linearized operator. Rather, 
we only say that it is only for the waves that do have the latter property that the numerical 
methods developed in this paper will be guaranteed to converge (with possible restrictions, 
as described in subsequent sections). For solitary waves whose corresponding linearized 
operators have more positive eigenvalues than the number of the wave's component, our 
modified CGMs are not guaranteed to converge. 

2 Comparative review of iterative methods for nonlinear prob- 
lem (11.41) and for linear problem (11.61) 

The reader who is not interested in such a review and wants to see the description of new 
algorithms proposed in this paper may skip to Section 3. 
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For reference purposes, let us rewrite (jl.4p as 

L(°)u = 



(2.1) 



An important role in our analysis will be played by the linearized operator L of the nonlinear 
Eq. (|2.ip . Let u be the exact solution of (|2,ip . u n be the approximation of that solution 
obtained at the nth iteration, and 

u n = u + u n , where \u n \ <C \u\. (2-2) 

Then 

L^u n = L^u + Lu n + 0{\u n \ 2 ) = Lu n + 0(j-u n | 2 ) . (2.3) 
For example, for the nonlinear equation (jl.5p . the linearized operator is 

L = -M + l/ (cos 2 x + cos 2 y) + 3n 2 . (2.4) 

In our discussion, the linear system (jl.6p is a counterpart of the nonlinear equation (|2.ip and 
matrix A is the counterpart of the linearized operator L. For Hamiltonian wave equations 
that give rise to (12. ip . operator L is always self-adjoint, which is the counterpart of matrix A 
being real and symmetric. Having stated this correspondence between the nonlinear problem 
(|2.1|) and its linear counterpart (|1.6j) . we will, from now on, refer to them interchangeably, 
because our statements will apply to both of them equally 

As we will explain below, the linearized forms of the iterative methods proposed in 
[J IS for solving the nonlinear problems (jl.ip and (jl.4|) are related to Richardson's 
method_| (see, e.g., [19], p. 274) of solving the linear problem (|1.6j) : 

Yn+i = y„ + (Ay n - b)Ar , (2.5) 

The necessary condition for the iteration scheme (|2.5p to converge to the solution, y = 
A _1 h, is that A be negative definite. Choosing the parameter At in (|2.5p so that At < 
— 2/A m i n , where A m i n is the minimum (i.e., most negative) eigenvalue of A, one ensures 
that Richardson's method for (|1.6j) with a negative definite A converges. 
A naive counterpart of (|2,5p for the nonlinear problem (|2.ip is 

u n+1 = u n + L^u n AT. (2.6) 

However, in most cases it will not converge to a solitary wave. Indeed, from (|2.3p . the 
linearized form of (|2.6p is 

u n +i = u n + Lu n Ar . (2.7) 

This equation is a counterpart of (|2,5p . and so a necessary condition for its convergence 
is that operator L be negative definite. However, for most fundamental single-component 
1 also known as Picard's or fixed-point iterative method 
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solitary waves, the corresponding linearized operator L has one positive eigenvalue (see the 
last paragraph of the Introduction), which causes divergence of small deviations u n in (12, 7p 
and, thereby, divergence of the iterations (|2.6p . 

To overcome this divergence, the iterative methods of [TJ [8] replaced L^°'u n in (|2.6p 
by a modified expression 

L (o,mod) Un guch that . ^ ( L (o) u = o) 44> (L(°< mod )u = 0) and 
(ii) the corresponding linearized operator L( mod ) was guaranteed to be negative definite. 
For example, in [8], 

L (0,mod) u= _ LL (0) U; (2.8) 

so that L( mod ) = —I? . In terms of the linear system (jl.6)l where ^4 is real symmetric but 
not sign definite, this is equivalent to applying Richardson's method (|2.5p to the so called 
normal equation 

- A 2 y = -Ab, (2.9) 

where the matrix on the left-hand side is always negative definite. 

In [7] , we proposed an alternative idea of constructing a neg ative definite L( mod ). This 
idea will play an important role in the development of this paper, and here we will explain it 
as applied to the linear system (11. 6h , As noted above, for most single-component fundamen- 
tal solitary waves, the linearized operator L has only one positive eigenvalue. Accordingly, 
let A^ 1 ) be the only positive eigenvalue of matrix A in (jl.6p . with the corresponding eigen- 
vector being In the context of the nonlinear equation (12. the eigenvector of L 
corresponding to its only positive eigenvalue is related to the solution u of (|2.ip |7], as we 
will explain later. Hence it can be considered as approximately known once the iterative 
solution u n is sufficiently close to u, which we will always assume to be the case (see ()2.2p ). 
Now, instead of (]1.6p . consider an equivalent system 

I- 7 z« (z (1) ) T ) (Ay-b) =0, (2.10) 



where I is the identity matrix and 7 is any number greater than one. Note that the term 
z (!) (zW) 71 is the projection matrix onto the unstable direction zW. Using the spectral 
decomposition of a real symmetric p x p matrix A: 

A = f>«z«(z«) T , ( 2.ii) 

i=i 

where the eigenvectors form an orthonormal set: 

z«) r z^ = { " . (2.12) 




it is straightforward to show that (|2.1(jp and hence (|1.6p is equivalent to 

^( mod V - [I - 7 zW (z (1) ) T ) b = 0, (2.13a) 



6 



where 

A (mod)^ A (i) (1 _ 7)z (i)^ z (i)j^ + £ A (O z (0( z «y . (2.13b) 

since A^(l - 7) < and A® < for t > 2 by our assumptions about A and 7, the 
matrix A( mod ) in dugED is negative definite, and hence Richardson's method (|2.5|) for it 
can converge. 

Note that if one chooses [7J 

7 = l + l/(A (1) Ar) (2.14) 

and uses the resulting ^4.( mod ) in Richardson's method (|2.5p . one thereby forces the projection 
of the error y n = y n — y on the direction of the unstable eigenvector z^ 1 ) to be zero at 
every iteration. This idea of zeroing out the projection of the iteration error on the unstable 
direction lies behind the original Petviashvili method [2] and its generalization for Eq. (jl.4p 
[7]. This is also the idea that we will use later in this paper for the development of modified 
CGMs. 

So far we have discussed how Richardson's method for (jl.6p with a sign indefinite A can 
be made to converge. The next issue is, how fast it converges. This is quantified by the 
convergence factor R, such that the number of iterations required for the iteration error to 
decrease by a factor of e is (asymptotically) inversely proportional to logi?. For the optimal 
value of At, the convergence factor of Richardson's method is (see, e.g., Sec. 5.2.2 in [20] , 
or [8j): 

_ cond(^( mod )) + 1 
R ~ cond(A( mod )) - 1 ' (2 " 15) 

where the condition number is related to the eigenvalues of A^ mod ^: 

cond(^( mod )) = A min /A max (2.16) 

(recall that all the eigenvalues of J 4( mod ) are negative and so cond(A( mod )) > 1). For 
cond(A) 3> 1, convergence to a prescribed accuracy occurs in 0(l/log(i?)) = 0(cond(74)) 
iterations. That is, the greater the condition number, the slower the convergence of the 
iterative method. A common method to reduce the condition number is to use a precondi- 
tioner, which can drastically lower |A m i n | (see, e.g., Lecture 40 in [19]). The preconditioned 
Richardson's method is 

y n +i =y n + B- 1 (Ay n -b)AT, (2.17) 

where B is the preconditioning matrbj^|. For the Richardson-type method (|2.6p with 
replaced by i(°' mod ) ) a convenient preconditioning operator, N, has a form similar to that 
of M. For example, if M is as in (|1.3p . the convenient form for N is: 

N = c-V 2 , O0, (2.18) 



2 For example, the well-known Jacobi and Gauss-Seidel iterative methods reduce to (|2.17|l for appropriate 
choices of B. 
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where c is some constant. 

Even if the condition number is lowered by reducing |A m i n | via preconditioning, the 
convergence of an iterative method may still be slow due to |A max | being close to zero. 
In the context of the nonlinear problem (jl.4p . examplified by (jl.5p . this occurs when the 
propagation constant \x is close to the edge of the spectral bandgap of the linear potential 
(i.e., Vo(cos 2 x+cos 2 y) in (jl.5p ); see, e.g., Fig. 4(a) in [21] or Fig. Q] in Section 6 below. In [9j 
[8] we proposed to accelerate the iterations by eliminating the slowest-decaying eigenmode, 
i.e., the component of the error u n that corresponds to the eigenvalue A max - This effectively 
removes this slowest eigenmode from the spectrum of L^ mod ^ and thereby increases the 
magnitude of its maximum eigenvalue. This, in its turn, increases the convergence rate of 
the method via (|2.16p and (|2. 15j) . The slowest mode can be eliminated in exactly the same 
way as the unstable mode (i.e., the mode with the positive eigenvalue AW) in (|2.1U|) . The 
slowest mode can be approximated by [9] 

Z SlOW K ^ _ ^ 

and a preconditioned Richardson's method is then applied not to (jl.6p but to 

/ - 7 zW (z (1) ) T - 7 slow z slow (z slow ) T ^ {Ay - b) = 0, (2.20) 

where 7 slow is computed similarly to (|2.14p . The explicit form of a counterpart of (|2.20p for 
the nonlinear problem (|2.ip will be given in Section 6. 

In [7J [8] we demonstrated that when the convergence of an iterative method is slow, the 
slowest mode elimination technique described above can considerably (by a factor of five to 
ten times) accelerate the convergence. However, it is well known (see, e.g., |19j . p. 341 and 
|20j . p. 451) that the fastest generic method for a linear system (|1.6p with a real symmetric 
and sign definite matrix A is the Conjugate Gradient method (CGM), whose algorithm is 
given by Eqs. (13. ip of the next section. The convergence factor of the CGM satisfies (see, 
e.g., [19], p. 299): 

R~ Vgjgg+l, (2 . 21) 
^cond(A) - 1 

which for matrices with large condition numbers implies considerably faster convergence 
than does (I2.15p . Therefore, if one could modify the CGM so that it would be guaranteed 
to converge when A (or the linearized operator L) is not negative definite but has one 
positive eigenvalue, then such a modified CGM is expected to be faster than a Richardson- 
type method accelerated by mode elimination. 

In this paper we present several versions of such a modified CGM for nonlinear problem 
(|2.ip and its generalizations and compare them with the Richardson-type methods of [7, 13] 
accelerated by mode elimination. We find that, as expected, the modified CGMs are the 
faster of these two groups of methods. As we noted after Eq. (12. 14ft , the main idea behind 
all these versions of the modified CGM is to eliminate the component of the error u n that 
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is aligned along certain unstable eigenmode(s) of the linearized iteration operator. Our 
modified CGMs are guaranteed to converge in almost all situations where a fundamental 
solitary wave with a prescribed (set of) propagation constant(s) is sought. Moreover, even 
for those rare situations where our guarantee is formally void, we will show a simple and 
practical way to still make the corresponding modified CGM converge. For the situations 
where a solitary wave with a prescribed (set of) power(s) is sought, our modified CGM is 
guaranteed to converge in the same parameter space where the Richardson-type method 
[131 122] converges. 

3 Modified CGM for solitary waves with prescribed propa- 
gation constant 

In this section, we will first state the algorithm of the CGM for the linear Eq. fjl .6[> . Then, 
we will briefly review the generalized Petviashvili method [7] for solitary waves and re-cast it 
in a form more convenient for the present analysis. Finally, we will propose a modification of 
the CGM for finding fundamental solitary waves with a prescribed value of the propagation 
constant. All analysis in this section is done for the single-component case; its generalization 
for multi-component solitary waves is presented in Section 5. 

3.1 CGM for the linear Eq. (ITBjl 

The original CGM algorithm for Eq. (jl.6p with a real symmetric and sign definite matrix 
A and starting with an initial guess yo, is: 

r = Ay - b, d = r , (3.1a) 

{ r ni d n ) 

a - = -Jld^-y (3 - lb) 

Yn+i = Yn + a n d n , (3.1c) 

r n+ i = Ay n+1 - b, (3. Id) 

= - <%■■■ ^ (3.1e) 

d n+ i = r n+ i + (3 n d n . (3. If) 

Here (a, b) denotes the natural inner product between the real- valued vectors a and b. 
Equations (|3.1ap define the initial residual vector Tq and the initial search direction do- 
Equation (|3.1cj) updates the solution by adding to the previous solution a vector a n d n along 
the search direction d n . The value of a n is set by ()3.1b|) to guarantee the orthogonality of 
the new residual, defined by (|3.1dj) . to the search direction d n : 

(r n+1 ,d n ) = 0. (3.2) 
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Finally, Eq. (|3.1fj) updates the search direction, where f3 n is set by (|3.1ep so that 

(d n+ i, Ad n ) = 0. (3.3) 

The condition that A be real symmetric is inherently required for the derivation of 
algorithm f)3. 1 [) . In particular, it guarantees that the orthogonality relations (|3,2p and (|3,3p 
between quantities at two consecutive iterations imply more general orthogonality relations 
among such quantities at all iterations (see, e.g., [23J or [19J, p. 295): 

(r n+1 , d j )=0, j<n, (3.4) 

(d n+1 , Adj) = 0, j<n. (3.5) 

The orthogonality relation (|3.5p with a sign definite A ensures that the CGM at each 
iteration produces a search direction that is linearly independent from all previous search 
directions. This implies that the error (in a certain norm) decreases with each iteration; this 
fact can be restated by saying that the CGM is an optimal iterative method for Eq. (jl.6p 
with a sign definite matrix A (see, e.g., |19j . p. 296). 

The condition that A be sign definite is also needed to guarantee that the denominators 
of a n and (3 n never vanish (or, for practical purposes, never become too close to zero). Thus, 
applying the CGM to an equation with a sign indefinite matrix A would come without the 
guarantee that this method would be optimal and that it would even converge. However, it 
should be emphasized that the CGM may converge when applied to problem (|1.6jl with a sign 
indefinite matrix A (see, e.g., [23], or [19J, p. 301). This is in stark contrast to Richardson's 
method (|2.5p . which, for a generic initial condition yo, is guaranteed to diverge in such a 
case. 

The CGM (|3.ip is straightforwardly extended to solve a nonlinear equation 

K {0) v = (3.6) 

whose linearized operator K is self-adjoint, by replacing Ay — b with K^'v in (|3.1aj) and 
(|3.1dp and A with K in (|3.1bp and (|3. lej) . In fact, for the CGM applied to a nonlinear 
problem, several choices of parameter (3 n that are equivalent in the linear approximation, 
will produce different versions of the method. Two most widely used such versions are 
known as Fletcher-Reeves and Polak-Ribiere (see, e.g., [M])- in the following we will not 
attempt to compare these two versions because the focus of this paper is not on the extension 
of the linear algorithm (13. ip to the nonlinear case but on the extension of that algorithm 
to the case when the counterpart of matrix A has one eigenvalue of the opposite sign than 
the rest of the spectrum. 

Let us also note that, unlike the algorithm (13. ip for linear equations, its extension to 
nonlinear problems described in the previous paragraph may fail if the linearized operator 
K is negative definite but has one or more eigenvalues close to zero. This can occur if, at 
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some iteration, the search direction d n becomes aligned primarily along the eigenmode(s) 
corresponding to the small eigenvalue(s). Then the denominator in (I3.1bj) will be small and 
hence the increment a n d n of the solution in (j3. lcj) . large. If the nonlinear equation admits 
more than one solution for the considered values of its parameters, which typically occurs 
near a bifurcation, then a large shift, a n d n , of the solution may cause the iterations to 
converge to a different solution than initially intended, or to diverge. Such a failure of the 
method can be avoided simply by changing the initial condition yo- We will discuss this in 
more detail in Section 6. 

3.2 Review of the generalized Petviashvili method 

In [7j, we proposed the following iterative method for finding the fundamental solitary wave 
u of the nonlinear problem f)2 . 1 p : 



where 



rt(\- 7 ^l, (3.7) 
(u n ,Nu n ) j 



(/, g) ^ J /(x)^(x)dx 

and the integration is over the entire spatial domain. The role of the 7-term is discussed at 
the end of this subsection. In (13. 7|) . a self-adjoint and positive definite differential operator 
N with constant coefficients is chosen so that u is an approximate eigenvalue of N^ 1 L 
corresponding to its largest eigenvalue A^: 

N^Lurs X W u. (3.8) 

The Sylvester law of inertia (see, e.g., [25]) guarantees that the signs of the eigenvalues of 
N^ 1 L are the same as the signs of the respective eigenvalues of L. Then, according to our 
definition of a fundamental solitary wave at the end of Introduction, A^ 1 ) may be positive 
and all the other eigenvalues of N _1 L are negative. 

The explicit form of ./V is usually taken to be similar to that of M. For example, 
when M is given by (|1.3p . N has the form (I2.18|) where the constant c can be computed 
algorithmically from u n [7]. (For the equation (11. ip with power-law nonlinearity, N = M 
and the equality in (|3.8p is exact [2].) The constant 7 in (|3.7p is given by (|2.14p and 
A^ 1 ) = An is found from 

A n x) = (u n , Lu n )/(u n , Nu n ) . (3.9) 

In practice, N, 7, and A^ 1 ) are computed until the iteration error reaches some predefined 
tolerance, after which their last-computed values are used for all subsequent iterations. Since 
N is a differential operator with constant coefficients, it has a simple representation in the 
Fourier space. Therefore, quantities like N~ l L^u n , Nu n , etc. are easily computed using 
the direct and inverse Fast Fourier Transforms, which are available as built-in commands 
in all major computing software. 
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In what follows we will frequently refer to the linearized form of (|3.7p . which is: 

u n+1 -u n = (N^Liin - 7 fej ^ t?j At . (3.10) 

Note that the operator N~ l L in the leading term of (|3.10j) is not self-adjoint. Therefore, 
we introduce a commonly employed change of variables to rewrite that equation in a form 
involving a self-adjoint operator. Namely, we denote 

v = N 1/2 u, v n = N 1/2 u n , K = N~ 1/2 LN- 1/2 , (3.11a) 
K (0) Vn = Ar-i/2 L (0) nn _ (3 . llb) 

Note that operator K is self-adjoint, since L is self-adjoint and N is self-adjoint and 
positive definite. Let us stress that K and v have been introduced only for the convenience 
of carrying out the subsequent analysis with self-adjoint operators. Once the derivation 
of the algorithm in this framework is complete, we will recast it in terms of the original 
variable u n and operator L. 

With the notations (|3"?TTj) . Eqs. (pTH]l . (|3"7T|) . and ([3TTH are rewritten as: 

Kv^X (1) v, (3.12) 

Vn + 1 ~Vn= - J A A T , (3.13) 

v n+ i -v n = (Kv n - 7 v) At . (3.14) 

Equations (|3.13p and (|3.14p are the nonlinear and linearized counterparts of the linear 
problem (|2.10p with a symmetric matrix A. The 7-terms in these equations are needed 
to eliminate the component of the error v n +i along the unstable "direction" v; see (|3.12p 
and the text after (|3.8p . If the equality in (|3.12p (and (|3.8p ) is exact, which occurs only 
for equations with power-law nonlinearity, (jl.ip . this elimination is complete. However, 
for the majority of equations (jl.4p with a more complicated nonlinear term, v is not an 
exact eigenmode of K and hence the 7-terms in (|3. 13[) and (|3.14p strongly suppress, but do 
not completely eliminate, the component along the unstable eigenmode. Nonetheless, this 
strong suppression is sufficient to turn the unstable eigenmode into a stable one. 

In the remainder of this section, we will continue using Eq. (|3. 14[) with the self-adjoint 
operator K for all derivations, and only at the final steps of those derivations will convert 
to the original variable u n . 

3.3 A modified CGM 

As we noted at the end of Sec. 3.1, the straightforward nonlinear generalization of the CGM 
applied to Eq. (|3.6p whose linearized operator K is self-adjoint but not sign definite, may 
diverge. For equations that we consider in Section 6 below, this method indeed diverges. To 
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modify the method so that it could converge, one can use the idea stated before Eq. (|2.8[) . 
Namely, replace (13.6P with an equivalent equation _ff(°< mod )t; = such that the correspond- 
ing linearized operator i^( mod ) is self-adjoint and sign definite. In Appendix 1 we show that 
the expression in (|3.15ap below, which mimics that on the right-hand side of (|3.13p . comes 
close to satisfying these properties: (i) it is equivalent to K^v, (ii) its linearized operator 
_?r(mod) is approximately self-adjoint, and (iii) its eigenvalues that are not too close to zero 
are all negative. Statements (ii) and (iii) are approximatqj rather than exact because so is 
relation (|3,12p . We will discuss the implications of the approximate character of statement 
(ii) after we present the algorithm of the modified CGM. 
Thus, denoting 

K (0,mod) w = K (0) v _ r (V, K(°)V) v 

\V, V) 

r = l + ^y, (3.15b) 

^ (mod) - = R ~ _ r (V, Kv) v 

(v,v) 

we propose the following modified CGM for finding the fundamental solitary wave of (|3.6[) . 
It mimics algorithm (|3.ip with the following modifications: Ay — b in (|3.1ap and (|3.1dp is 
replaced with K (°> mod )v, A in (13. lb[) and (]3. lef) is replaced with K^ mod \ and vectors r and 
d are renamed r and d. Upon the change of variables 

f n = iV 1/2 r„, d n = N l / 2 d ni (3.16) 

and ([5XT]) . this modified CGM is: 

r = N^L^^uo, d = r , (3.17a) 
(Nr n , d n ) ,„ 17U 

= - W-Od.,*.) - (3 ' 17b> 

u n+ \ =u n + a n d n , (3.17c) 

r n+1 = N-W^Un+i (3.17d) 

(r n+1 , L^d n ) 
Pn ~ (L(™*)d n , d n ) ' ( '- 17ej 

d n +i = r n+ i + (3 n d n . (3.17f) 



where 



L (0,mod) u = L (0) u _ p {U,L(°) U ) ^ 

(u, Nu) 

L ^od) d = Ld _ T (u, Ld) 

(u, Nu) v ' 



3 However, as we point out in Appendix 1, in all the examples that we tried and some of which are 
reported in Section 6, statement (iii) holds in a more definite sense, i.e.: all eigenvalues of K^ ^ are 
negative. 
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In a practical implementation of algorithm (|3.17p . iterations should be carried out with 
the generalized Petviashvili method (|3.7p until the error reaches some predefined tolerance 
(usually, between 1 and 10%), so that the computed parameters of operator N and the 
constants AW and r may be considered as known with the corresponding accuracy. After 
that, the iterations should be carried out with the modified CGM (|3.17p . which is expected 
to converge considerably faster than the generalized Petviashvili method. 

Operator x( mod ) in ()3. 15cj) . and hence the quadratic form (d n , K^ mod 'd n ) in the coun- 
terparts of (|3.17b .e). would be guaranteed to be negative definite if in relation (|3. 12j) (or, 
equivalently, in (|3.8p ). the equality held exactly. However, this takes place only for equa- 
tions with power-law nonlinearity, (jl.ip . Therefore, for some nonlinear wave equations, 
(dn, K( mod )d n ) could be sign indefinite and, importantly, arbitrarily close to zero. In such 
situations algorithm (|3.17j) could fail, and we will now point out when this should be ex- 
pected. As we explain in Appendix 1, the slight non-self-adjointness of if( mod ) can cause 
the quadratic form to become too close to zero when some of the negative eigenvalues of 
the linearized operator L of (12. ip are close to zero. In single-component equations, this typ- 
ically occurs when the solitary wave bifurcates from the edge of the continuous spectrum, 
while for multi-component equations this can also happen when, say, an asymmetric wave 
bifurcates from one with a symmetry (see, e.g., [26] and references therein). 

Our numerical experiments confirm that algorithm (13.17P converges whenever the eigen- 
values of L are sufficiently far from zerc@, and also that is may diverge in the opposite sit- 
uation, as described above. Recall, from the last paragraph of Section 3.1, that in exactly 
the same situation, a CGM for a nonlinear problem can fail even if the quadratic form 
(d n , K^ mod ^d n ) were negative definite. Thus, there are two different reasons that can cause 
the modified CGM to fail, but, fortunately for the applications of the method, they both 
can do so in the same situation, namely, when L has near-zero eigenvalues. Moreover, they 
both require that at some iteration, the search direction become primarily "aligned" with 
the eigenmode of L corresponding to a small eigenvalue. Therefore, a simple and practical 
way to avoid such divergence is to change the initial condition in a certain manner. In 
Section 6 we show that this way does indeed work. 

Before concluding this section, let us note that one could develop a modified version of 
the CGM based on a different idea than the elimination of the unstable eigenmode from 
the underlying operator, as in f)3. 15|) . Instead, one can employ the original operators 
and K but force all search directions d n to be orthogonal to the unstable eigenmode of K 
approximated by the exact solitary wave v. Since (|3.12p holds approximately, the search 

4 Note that if the solitary wave is translationally invariant, L has a zero eigenvalue corresponding to 
the respective translational eigenmode. As we show in Appendix 1, such a mode presents no problem for 
convergence of the iterations since it would simply slightly shift the solitary wave. Therefore, in the following 
we ignore the possible presence of such zero eigenvalues of L. 
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directions can be made only approximately orthogonal to the true unstable eigenmode of 
K. One can show that this causes this other modified CGM to be prone to failure under the 
same circumstances as the modified CGM ()3.17|) . i.e. when L has eigenvalues that are close 
to zero. Our numerical experiments show that the simple trick that can help algorithm 
(|3.17p overcome that failure, is not as effective for the other modified CGM. Moreover, the 
coding of the latter method is slightly more complicated than that of (|3.17p . For these 
reasons, we only present that alternative modified CGM in Appendix 2, but do not discuss 
it further in this paper. 



4 Modified CGM for solitary waves with prescribed power 

We first review the preconditioned imaginary-time evolution method (ITEM) [13], which 
is a Richardson-type iterative method for finding solitary waves with a specified value of 
power, and then present a modified CGM that converges under the same conditions as the 
ITEM, but faster. 

4.1 Review of the ITEM 

When one seeks a solitary wave with a prescribed value of power (II. 7ft . problem (12. 1|) 
should be redefined because the propagation constant /i, which enters via operator M 
(see (|1.3|) ). is no longer known exactly. Instead, it is estimated at each iteration, as shown 
below. Thus, we rewrite (|1.4p as 

L^u-/j,u= 0, (4.1a) 

where 

L(°% = V 2 n + F(n,x), „ = {fi ^ f ^ , (4.1b) 

(/(«)>«) 

and f(u) is any function of u; its choice will be specified shortly. 

The preconditioned ITEM whose convergence conditions are found in [13] is: 

^ = -^Tl ' 4 - 2a 

(N 1 u,u) 

Un+l = u n + N- 1 (l^U n - fl n U n ^j At , (4.2b) 

u n+ i = u n+1 A— — -. (4.2c) 

Y {u n +l,U n+ i) 

Here P is the prescribed value of the power and N is a preconditioning operator of the 
form (|2.18p where now c is an arbitrary (unlike in Section 3) positive number. 

In what follows we will need a few facts about the linearized method (|4.2p |13j . First, 
the condition that (u n ,u n ) = P = const entails the orthogonality between the error and 
the exact solution: 

(u n ,u) = 0{u 2 n ) . (4.3) 
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Next, the linearized form of the operator in parentheses in (|4.2b[) is: 

Lu n = Lu n 7T— i r-u . (4.4) 

{N~ 1 u,u} 

(Note that operator L involves the propagation constant [i. Even though \jl is not specified 
when we seek the solitary wave, we can still use it in analyses of iterative methods.) Using 
the last two equations, one can straightforwardly show that the last line, (14.2cp . of the 
ITEM does not change the linearization of the preceding line. Equation (|4.2cj) is needed 
to ensure that the power of the solitary wave equals P exactly rather than in the linear 
approximation . 

The ITEM (|4.2|) converges when its linearized operator, has only negative eigenvalues 
(see the footnote at the end of Section 3). Loosely speaking (see [13] for a more precise 
statement), for a wide subclass of equations (|1.4p . this occurs when operator L has none or 
one positive eigenvalues. Moreover, in the latter case, the following condition has to hold: 

dP/dfi > . (4.5) 

Below we will use the change of variables (|3.11a|) and similar notations: 

K m Vn = N -vi L m Un , K m VN = K m Vn _ ^ nN -\ n , (46) 

where \i n is given by ([4.2aj) . The linearization of K.(°>v n is an operator 



When the convergence condition of the ITEM (|4.2p . stated near (|4.5p . hold, operator K, is 
[13j self-adjoint and negative definite on the space of functions satisfying the exact form of 
the orthogonality relation (|4.3j) : 

(v n ,N- 1 v} = 0. (4.8) 

4.2 A modified CGM 

The idea of this modified CGM is to ensure that both the iteration error v n and the search 
direction be orthogonal, in the sense of (|4.8p . to the exact solution v. Then, according 
to the last sentence in the previous subsection, operator JC is guaranteed to be negative 
definite on this space of functions, and hence the following modified CGM is guaranteed to 
converge under the same conditions as the ITEM (|4.2|) : 



f = X: ( %o, d = r -^(N- 1 v ,r )v , (4.9a) 
a n = - / f n ^ n J, , (4.9b) 




Vn+l = V n + Ct n dn , V n +1 = V n +1 \ I 77 77 T7 T , (4.9c) 

w (V n +1,N i V n +l) 
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f n+1 = K^v n+l (4.9d) 

(f n+ i, Kd n ) - —{v n+ i,lCd n )(N~ 1 v n+ i,f n+1 ) 

Pn = P pjr , (4.9e) 

dn+i = r n +i + Pndn ~ p(iV _1 i; n+ i, f n+1 + p n d n ) v n+1 . (4.9f) 
The definitions of do in (|4.9ap and c2 n +i i n (|4.9f| ensure that 

(4,JNT 1 «n) = (4-10) 

at every iteration. Here we use the fact that our iterative solution has a fixed value of 
power: 

(v n ,N- 1 v n ) = P (4.11) 



for all n. Since the error v n at the nth iteration satisfies the orthogonality condition (|4.8p . 
the first of Eqs. (|4.9cj) and Eq. (|4.10p ensure that the error £ n +i at the next iteration also 
satisfies (14.8(1 . The second of Eqs. (14.9c)) does not change this relation; it only forces (|4.1ip 
to be satisfied exactly rather than in the linear approximation. Equations (I4.9b -d) entail 
the counterpart of (|3.2|) : 

<f n+ i, d n ) =0(v 3 n ), (4.12) 
and Eqs. (|4.9e .f) entail a counterpart of (13.3|) : 

JCd n ) = 0{vl) . (4.13) 

We demonstrate these relations in Appendix 1. 

Let us note that condition (|4,10p . in addition to ensuring (|4.8p at each iteration, playes 
one other important role in this algorithm. Namely, it ensures that in the inner product on 
the left-hand side of (I4.13p . /C is a self-adjoint operator: see the sentence before (I4.8p . This 
property of K, entails the counterpart of (|3.5p , which together with negative definiteness of 
K. implies that all search directions d n are linearly independent and the error decreases as 
the iterations proceed. In Section 3.1 we noted that this means that the CGM is an optimal 
iterative method (see, e.g., [19], p. 296, and |24J). 

Thus, the modified CGM (|4.9p is guaranteed to converge under the same conditions 
as its Richardson- type counterpart (14. 2p . This, mathematically, is a more rigorous result 
than what we obtained for the modified CGM (|3.17p . which is guaranteed to converge in a 
slightly narrower parameter space than its Richardson- type counterpart (|3.7p . 

Finally, we present algorithm (|4,9p in the original variables: 

r = N- 1 £<®Uo, d = r - ^(u ,r )u , (4.14a) 

-TOT ( " 4b> 
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u n +i =u n + a n d n , u n+ i = u n+ i , 



P 



(u n+1 ,u n+1 ) ' 



r n+1 = N^C^Un+i 

r n+ i, Cd n ) - —(u n+ i,Cd n )(u n+ i,r n+ i) 
(d n , Cd n ) 



1 



d n +i = r n+ \ + {3 n d n - —(u n+ i, r n+ i + f5 n d n ) u. 



P 



n+l 



Here C, which is the linearization of 



jr(o) =L (oo) u {N-W,L.m Un ) 
(N l u n ,u n ) 



is given in (|4.4p . 



(4.14c) 
(4.14d) 

(4.14e) 
(4.14f) 

(4.15) 



5 Modified CGMs for mult i- component solitary waves 



Here we present the multi-component versions of the modified CGMs ()3.17j) and (|4,14p for 
fundamental solitary waves with prescribed values of, respectively, propagation constants 
and quadratic conserved quantities that are the counterparts of power (jl.7p . We will not 
give derivations of these methods because they are similar to those found in Sections 3 and 
4. 

We begin with stating the generalized Petviashvili method from [7| for finding a multi- 
component solitary wave u = [u^\ . . . ,u^] T of the equation 



L(°)u = 



(5.1) 



whose linearized operator is L: 



u n +i — u n + 



(k) 



P (fe) T (°)n 



- — . / 0) TVT ( fc ) 

k=l \ ; J- ^ G n 



At. 



(5.2) 



Here N = diag (N^\ . . . , N^) where each has a form similar to the linear differential 
part of the kth component of L/ ); e.g., for a system of nonlinear Schrodinger-type equations, 



N (k) _ c (k) _ b (*)y2 j 



(5.3) 



and c^ k \ are computed by formulas found in [TJ. Furthermore, are the approximate 
eigenvectors of N _1 L: 

Le ( fc ) ra A^Ne^ , (5.4) 



which are mutually orthogonal with weight N: 

(eW.Ne^) = for k ^ m. 



(5.5) 
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Note that (|5.4p is the multi-component counterpart of (|3.8p . The form of e'^'s is related 
to the solution of (|5.ip by 

e«=u n; e« = (^V 1 ),...,p^V s )) T , (5.6) 

where pC 5 '*) = 1 and the rest of are found as explained in [7]. Finally, the constants 
and 7 (fc) in ([Ojl are found similarly to flEI]) and (I2TT4T) : 

A (fc) = ( e ( fc ),Le (fc) )/(e (fe) ,Ne (fc) ), 7 (fc) = 1 + l/(A (fe) At) . (5.7) 

The multi-component version of the modified CGM of Section 3 looks as (]3.17p . with the 
scalar functions u n , r n , d n being replaced by the S-dimensional vectors u n , r n , d n , operator 
iV being replaced by its matrix counterpart given before Eq. (|5.3p . and the operators L(°' mod ) 
and L( mod ) being given by expressions generalizing (|3.18p : 

s (e { n k \ L(°>u n \ 
L (o, mod ) Un = L (o) Un ^Ne« , (5.8a) 



fc=i 



r« = i + ^ , (5.8b) 



5 (ei fc) , Ld n 

L(-° d )d n = Ld„ - E r(fc) yV) — m Nei n ] ■ ( 5 - 8c ) 



k=i (e^ ; ,Ne^y 

A sample code of this algorithm for a two-component system considered in Section 6 is 
presented in Appendix 3. 

Let us note that the purpose of the terms under the sum in (]5.8p is to eliminate the 
eigenmodes e^^ of L whose eigenvalues could be positive. By our definition of the fun- 
damental solitary wave, stated at the end of Introduction, eliminating S such eigenmodes 
should be sufficient to turn the positive eigenvalues of L into negative eigenvalues of L( mod ). 
As explained in Section 3, due to the approximate character of relations (|5.4p . this unstable 
mode elimination still may not always guarantee convergence of the multi-component mod- 
ified CGM (|3.17p . (|5.8p when L has small negative eigenvalues. However, a simple approach 
based on shifting the initial condition, as discussed in Section 6, overcomes this divergence. 

We now turn to the case where a set of quadratic conserved quantities related to the 
solution components' powers are prescribed. We need to introduce a number of new no- 
tations. Let us denote the s-component vector of these quantities by Q, so that its kth 
component is 

Q(fe) = ^ 9 («) P (O j p(0 = („(0 >u (0) ? k = l,...,s<S, 1 = 1,... ,S. (5.9) 
i=i 

Note that the number of these conserved quantities can be less than the number of the 
components of the solitary wave: s < S. To emphasize this fact, we use a different vector 
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notation for Q than for u. The number of the components of Q equals the number of the 
propagation constants that one can prescribe independently. One example of this situation 
is the system of three waves interacting via quadratic nonlinearity [27] (see also [8] ) . Another 
example, which we discuss in more detail below to make our notations clear, is the system 
describing the evolution of pulses in a two-core nonlinear directional coupler 

m$ k) + ui k J + (V^| 2 + k\U^\ 2 ) u& + = 

.^(fc+2) + ^fc+2) + ^(* +2 )|2 + K \u(k)^ ^(fc+2) + ^(5-*) = 

where and [/( fc+2 ) are the orthogonal polarization components of the pulse in the kth 
core. Upon the substitution 

/ l/V>(x,t) \ ( uW(x) \ 



k = 1,2. 



(5.10) 



(2), 



u (3 Hx) 



V 



1/ 



(4), 



Ci)t 

C2) t 



(5.11) 



where can be chosen to be real-valued, system (|5.10|) reduces to: 



^ + ((«( 2 ») 2 + ^)>< 2 »+«w 

ni 3 j + ((u^) 2 + K (u«) 2 ) n( 3 ) + 





( 







\ 




/ 





\ 





























u (3) 




M = 











\ 







/ 




V 





/ 



(5.12) 



where /!=(// 



Thus, in this example, S = 4, s = 2, and the vector of quadratic 



conserved quantities (verified from the time-dependent equations (|5.10p ) is 

/ pd) \ 



Q 



p(l) + p(2) 
p(3) + p(4) 




p(2) 
p(3) 
p(4) 



(5.13) 



Generalizing the above example, one can see that the counterpart of Eq. (|5.ip when Q 
rather than fl is prescribed, is the following extension of Eq. (|4.1a|) : 

L(°°)u - U (N- 1 ^, U)- 1 (NTfy L(°°)u) = , 



(5.14a) 



Ou 



(5.14b) 

For example, in (|5Tl2"j) . i/ 00 ) u is the first term (the 4x1 vector), and £Y is the first factor 
of the second term (the 4x2 matrix) on the left-hand side. 

For the convenience of subsequent notations, we will assume that the matrix (q^) 
in (|5.9|) is put into the reduced echelon form (see any textbook on undergraduate linear 
algebra) and, in addition, its columns are rearranged so that 



Jkl) 




(5.15) 
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For example, for the 2x4 matrix in (|5. 13|) this would imply switching the second and third 
columns. Then for the powers of the solitary wave components whose indices equal the 
indices of the pivot columns of one has from (j 5 . 9 j) : 

S 

p(k) = Q (k) _ £ q {ki) p {i) ; k = l,... ,s<S. (5.16) 

l=k+l 

We will use this fact in the next paragraph. 

With the convention (|5.15p . the multi-component version of the ITEM is: 

u n+ i = u n + N- 1 (L(°°)u n - U n {N^Un, Un)- 1 (N- 1 ^,, L( % n }) At , (5.17a) 



(k) _ 4k) 
u n+l — u n+] 

where 



n(fc) _ v 5 n( kl ) P w 

' ^=^ q _^±i, k=l,...,s<S, (5.17b) 



n+l 



p(k) = ,Ak) , (fe) > , _ 1 „ 



In analogy to the algorithm (|4.2p for single-component equations, Eq. (|5. 17b|) does not 
change the linearized form of (I5.17ap ; its role is to guarantee that the s components of 
vector Q equal their prescribed values exactly rather than in the linear approximation. 

To our knowledge, the multi-component ITEM (|5.17p has not been presented in the liter- 
ature before; however, its close "relative" (related to a family of squared-operator methods) 
was stated in [8j. The convergence conditions for (|5.17p are [22]: (i) the s x s Jacobian 
dQ/djl must be nonsingular, and (ii) the number of positive eigenvalues of this Jacobian 
must equal the number of positive eigenvalues of the linearized operator L of Eq. (|5.ip . 
These conditions are a counterpart of (|4,5p . They guarantee that the linearized operator, 

Cu n = Lu n -U (N -1 W, U)- 1 (N -1 W, Lu n ) (5.18) 

of the expression 

£(0) 

u ra appearing inside the parentheses in (I5.17ap . is negative definite on 
the space of vector functions satisfying an analog of the orthogonality relation (|4.3p : 



(U, u n ) = 0. (5.19) 

The corresponding multi-component version of the modified CGM (|4.14p is: 

r = N-^^uo, d = r - U Q (Uq, U Q )~ l (U , r Q ) , (5.20a) 

(r n , Nd w ) 

a " = -(d^)' (5 - 20b) 



, (fc) *(fc) 

u n+ i = u n + a n d n , = 



n+l 



\ P { 
\ n 

(5.20c) 

r n+1 = N-^^V! (5.20d) 
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a (rn+i; £d ra ) - (Cd n , U n +i){U n +i, U n +i) 1 (Un+l, r n+ i) 

d n+ i = r n+ i + f3 n d n - U n +i (Un+i, W n+ i} -1 (Z4+1, r n+ i + (3 n d n ) . (5.20f) 

This method is guaranteed to converge to a fundamental solitary wave with a prescribed 
vector of conserved quantities Q under the conditions stated before Eq. f|5. 18jl . 



6 Numerical examples 

Here we compare the performance of the modified versions of the CGM presented in Sections 
3-5 with the performance of the corresponding Richardson-type methods accelerated by the 
mode elimination (ME) technique [9j [8] . 

The model equations are (jl,5p and its two-component extension: 



(6.1) 



V 2 u« + y (cos 2 x + cos 2 y)u« + (f^ 1 )) 2 + F^Hu^) 2 ) u« = 
VV 2 > + y (cos 2 x + cos 2 y)u^ + (f<W( u W) 2 + F^(u^) 2 ) = ^ 2 V 2 > 
where the constants 

fW = i, F^ = A, f( 12 ) = o.5. 

Note that Eq. (jl.5p is equivalent to Eqs. (2.1), (2.2) of [28], where A*[28] = 2Vo — /^this paper- 
The values of Vq and the propagation constant(s) or power(s) are specified below, so that 
each modified CGM and the corresponding Richardson- type method accelerated by ME are 
tested for three cases: mildly numerically stiff, stiffer, and stiffest. (In the stiffest case the 
methods take the most number of iterations to converge.) A code for the most complicated 
of these methods — the modified CGM (|3.17p . (|5.8p for the two-component Eqs. (|6.ip - 
is presented in Appendix 3. (A counterpart of this code for a single-component equation 
is considerably simpler and can be easily written based on a code from Appendix 1 in [7]. 
Codes seeking solitary waves with prescribed power are also much simpler than the code in 
Appendix 3 below, because they do not compute quantities N, e^ k \ etc.) 

We begin by describing the simulations setup for the methods where the propagation 
constant (s) is (are) specified. We compare the modified CGMs (I3.17P (for the single equation 
(|1.5p ) and (|3.17p . (|5.8p (for the two-component system (|6,ip ) with the ME- accelerated 
generalized Petviashvili method. The latter method for a single equation is [9]: 



u n +i = u n + 



TV-MO),, >n, Li°) Un ) lQW (C lQW , £ (0) H n) ..slow 
1 V J-J Un. T , , - v "»n. 7n , , „i , - , „i 



7 (u n , Nu n ) Un 7n N<flr) rn 

which for a multi-component system extends to: 



At ; (6.2a) 



u n +i — u n + 



4 fc) , L(°)u n \ /a slow T (°)n 

N -1 L (0) U _ V-yW-^U - -.slow V^n i^'^n/ 



slow 



At. 
(6.2b) 
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Here the notations are as in (|5.2p above and in addition, 

h /<j> s l° w T *slow\ 

^slow _ _ ~,slow _ -i _| IL \slow _ V^n > ,n o\ 

- u n 7n + A^ 1ow At' n ~ (*^ low , N*^ low ) ' 

where /i is the fraction of the slowest-decaying mode that is subtracted at each iteration. 
In [9] we advocated subtracting around 70% of such a mode for robust and nearly optimal 
performance; accordingly, here we use h = 0.7 in all examples. 

For all the methods listed in the previous paragraph, we designate the following sets of 
parameters in Eqs. (|1.5p and (|6.ip as corresponding to the mildly stiff, stiffer, and stiffest 
cases. 



For (fL5|) : mildly stiff 


^> 


^o = 4, 


^ = 


5.03; 






stiffer 


=> 


^o = 4, 


^ = 


4.95; 






stiffest 


=> 


^0 = 6, 


/" = 


7.89; 




For (HU) : mildly stiff 


V 


= 4, 


M « = 


5.03, 




= 5.5 


stiffer 




= 4, 




4.95, 




= 6.5 


stiffest 


=> v 


= 6, 




7.89, 




= 8.5 



(6.4a) 



(6.4b) 



Figure Q] shows schematically the power- versus- /i plots for Eq. (jl.5p . with the three cases 
of (|6.4ap labeled. Figure [2] shows the solution for the stiffer case for Eqs. (|6.1|) : note the 
vastly different amplitudes and widths of the two components. In the other two cases the 
solution looks qualitatively the same. The obtained solutions of the single Eq. (jl.5p look 
qualitatively as the first component of the solution in Fig. [2j In general, the closer the 
propagation constant to the edge of the band gap, the broader and lower the corresponding 
solution. 

For both the modified CGMs and the ME-accelerated generalized Petviashvili methods, 
we start the iterations by the non-accelerated generalized Petviashvili method, and when 
the error, defined as 

^ ((K% n )W, (LWu n )W) 

£n ~ h U {k) u {k) \ ' 

k=l \ u n , <Mi I 

reaches a threshold 



£ accele r ationt hreS hold = 5-10- 2 , (6.6) 

we begin the acceleration by either the modified CGM or ME and carry on the itera- 
tions until the error reaches 10~ 10 . The parameters of N, e^ k \ X^ k \ and 7^ (or their 
single-component counterparts) stop being computed at the threshold (|6.6p and the latest 
computed values of these parameters are used from that moment on. (As was demonstrated 
in [7j, for Eqs. (|6.ip the eigenvector e^ 2 ) of N -1 L is computed not too accurately. Yet, we 
show below that this does not prevent the modified CGM (|3.17p . (|5.8p from performing 
quite well.) 

In regards to the specific numeric value on the right-hand side of (|6,6p . we observed 
that the accelerated methods converge in all cases when this value is not too high (say, is 
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less than 1CT 1 ). If the error is substantially higher, the computed parameters of N etc. 
may be too inaccurate, so that the corresponding operator K( 0,mod ) would be "too far" 
from self-adjoint, which might then prevent the convergence of algorithms f)3. 1T|) or (|5.8p . 
In addition, only for the modified CGM in the stiffest case, the value of the acceleration 
threshold should not be too low (e.g., should be higher than about 5 • 10~ 3 ); otherwise, 
the modified CGM in this case would diverge. The explanation of this fact was sketched 
in the last paragraph of Section 3.1. Here we present specific details for it; for the sake 
of convenience, we do it for the single Eq. (|1.5p and use the transformed variables (|3.1ip . 

(pngp . 

First, as is well known (see, e.g., |28] and references therein), when ju in (jl.5p is very close 
to the edge of the spectral gap of the linear part of K^> (i.e. in the "stiffest" case considered 
here), the shape of the exact solution of (|1.5|) is very similar to that of the eigenmode(s) of 
K with negative near-zero eigenvalue(s). In particular, the solution is broad, as illustrated 
in the top panel of Fig. [2] Next, from (|3.17a .b) and an analog of (|2.3|) it follows that at the 
first CGM iteration, the step size «o along the search direction do is 

(K^ mod "iv , K^ mod ">v ) _ (if( mod )v , K (mo ^v ) 
a ° ~ (jf( mod ) if (°> mod )?; , K(°> mod ho) ~ ((if( mod )) 2 vq, jf( mod ) vq) > 

(Here the iteration's number is referenced from the start of the CGM algorithm.) The last 
strong inequality follows from two facts: (i) operator K (= N~ l l 2 LN~ 1 / 2 ) has a band of 
negative eigenvalues that comes very close to zero (see, e.g., Fig. 4(a) in [21J for a related 
problem), and (ii) it is precisely a superposition of the corresponding eigenmodes that 
dominates the error vq obtained by the generalized Petviashvili method after sufficiently 
many iterations (because these modes decay the slowest in Richardson- type methods; see, 
e.g., |8j). Having a?o 3> 1 and hence a large first step in the CGM adds to vq some superpo- 
sition of near-zero eigenmodes of K and thereby makes the solution at the next iteration, 
x>\, more broad and hence even more resembling the near-zero eigenmodes. Then the step 
size ax, computed at the next iteration, becomes even greater than cto, and the solution at 
the following iteration, V2, is made even more broad. This quickly leads to divergence of 
the iterations. (Practically, the solution converges to a nonlocalized two-dimensional plane 
wave.) On the other hand, if the acceleration by the modified CGM starts when the error 
vq is not too small and hence still contains a significant contribution from the eigenmodes 
whose eigenvalues are not too close to zero, the corresponding step size ao is not too large, 
and subsequent CGM iterations are able to gradually reduce the error v n . 

The above consideration explains why the CGM iterations should start when the itera- 
tion error is not too low. In practice, however, this is not a limitation, since in numerically 
stiff cases, one wants to start the acceleration by the CGM sooner rather than later since 
this considerably reduces the computational time, as we will demonstrate below. 
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The initial condition in all cases that we consider is taken as 

u = 1.5 • e~ {x2+y2) • (1 + e x x + e y y), e x = 0.1, e y = -0.2 (6.8a) 
for Eq. (jl.5p and as 

n c \ 

. e -(* 2 +?/ 2 ) . (l + €x x + eyy), e x = 0.1, £^ = -0.2 (6.8b) 





o 

for Eqs. (|6.ip . The exact values of the amplitude(s) and width of the initial condition 
are not essential; all methods converge for a wide range of these parameters. As for the 
asymmetry parameters e x and e y , these can be also varied quite a bit as long as the shape 
if the initial condition resembles a pulse with one main peak. However, these parameters 
should not be set to zero (or too close to zero) for the modified CGMs (|3.17|) and (|5.8p in 
the stiffest case. (That is, in the other two, less stiff, cases, the modified CGMs converge 
even for a symmetric initial condition.) The reason for this is similar to the reason why the 
acceleration threshold in (|6,6p should not be chosen too small. Namely, iterations of the 
generalized Petviashvili method which start with a symmetric initial condition, i.e., (16. 7p 
with e x = e y = 0, produce symmetric solutions at all subsequent iterations. These solutions 
increasingly resemble the near-zero eigenmode of the linearized operator L. Then when the 
CGM starts, the first step «o along the search direction do becomes too large, which leads 
to divergence of the iterations, as explained after Eq. (|6.7p . On the other hand, having 
an asymmetric initial condition leads to the iteration error being asymmetric, and such an 
asymmetric error can have a considerable content of eigenmodes of L whose eigenvalues 
are not close to zero. This reduces the step size ao at the first CGM iteration, and the 
modified CGM is able to converge. Note that from the practical perspective, starting with 
an asymmetric initial condition is not a limitation of the method. 

While the ME-based acceleration does not require an asymmetric initial condition for 
convergence, we still use the same expressions (I6.8P for all methods, so as to make our 
performance comparison uniform. 

We now comment on the choices of the "fictitious time" step size At. For each sim- 
ulation, we first selected its nearly optimal value for the corresponding non-accelerated 
generalized Petviashvili method. This takes just a couple of trials since the optimal At 
is only slightly less than the maximum value of this parameter for which the generalized 
Petviashvili method still converges. (See, e.g., Eqs. (2.5) and (2.6) in [9] and Fig. 1(d) in 
[8], based on the same equations, although for another iterative method.) Once this Ar opt 
has been determined, we use a slightly smaller (specifically, by 0.1) value of the step size 
for the method accelerated by ME; the analysis of ME [9] predicts that this should lead 
to a more robust and smooth peformance than using Ar op t. For example, if we find that 
Ar opt = 0.9 for the non-accelerated generalized Petviashvili method, we then use Ar = 0.8 
for the ME-accelerated version of this method. For the method accelerated by a modified 
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Iterative method 


Case in (|6.4a|) 


non-accelerated (|3,7p 


accel. by ME (|6.2a|) 


accel. by COM (|3.17l) 


mildly stiff, At = 1.1 


300 iterations, 130 s 


110 iterations, 60 s 


60 iteratons, 40 s 


stiffer, At = 1.1 


920 iterations, 390 s 


290 iterations, 170 s 


100 iteratons, 60 s 


stiffest, At = 1.0 


3700 iterations, 1560 s 


430 iterations, 240 s 


200 iteratons, 120 s 


Table 1: Comparative performance of the generalized Petviashvili method accelerated by 
ME (|6.2a|) and the modified CCM (I3.17|) for the single-component Eq. (fL5|). 




Iterative method 


Case in (l6.4bj) 


non-accelerated (|5.2p 


accel. by ME (|6.2bl) 


accel. by CGM 


mildly stiff, At = 1.0 


330 iterations, 310 s 


120 iterations, 130 s 


70 iteratons, 70 s 


stiffer, At = 1.0 


780 iterations, 740 s 


200 iterations, 220 s 


130 iteratons, 130 s 


stiffest, At = 0.9 


3330 iterations, 3110 s 


550 iterations, 610 s 


240 iteratons, 250 s 



Table 2: Comparative performance of the generalized Petviashvili method accelerated by 
ME (^2bl) and the modified CGM (|3"T7j) . (ffTKD for the two-component Eq. ffTTD . 



CGM, we also start the iterations using At = Ar op t — 0.1, although this does not noticeably 
affect its performance (since the CGM itself does not involve At). 

Finally, the domain for all our numerical simulations is [— 67T, 67t] x [— 6tt, 6tt], with 
2 8 x 2 8 grid points. 

Tables 1 and 2 list the numbers of iterations and time (both rounded to the nearest ten) 
for each of the three methods (non-accelerated generalized Petviashvili method and its two 
versions accelerated by ME and CGM) in the three cases (|6.4p for Eqs. (jl.5p and (|6.ip . In 
accordance to the note above, the value of At is listed only for the non-accelerated method. 
In Fig. [3] we plot the iteration error (|6.5p versus the iteration number for the stiffest case for 
Eq. (jl.5p . The error evolutions shown there are representative of all the other cases, except 
that in the less stiff cases, the curves corresponding to the ME and CGM are less jagged. 

We also verified that a straightforward extension of the CGM, as described at the end of 
Section 3.1, would diverge for all cases listed above and even for non-stiff cases (not listed 
here) . 

From these tables we see that, as expected, both ME- and CGM-based accelerations 
considerably improve the performance of the iterative method, with the stiffer the case, the 
more the improvement. Furthermore, the CGM-based acceleration performs better than 
the ME-based one by a factor of 2-2.5; again, the stiffer the case is, the more improvement 
the CGM provides over the ME. 

We now present similar results when the power (jl.7p . or its multi-component counterpart 
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(|5.9p . of the solution is specified. The non-accelerated methods are the ITEMs (|4.2j) and 
(|5.17p applied to the single- and multi-component Eqs. (|1.5j) and (|6.ip . respectively. When 
the ME-based acceleration is applied, the lines updating the intermediate solutions u n +i 
and u n+ i of these methods become: 



U r , 



N~ 



7 



slow 



(<^ low , N<$> w ) 



A slow 



At: 



(6.9a) 



u n +i 



u n + 



w , L(°)u n > 



slow 
In 



<j>slow 



^slow 

, N*^ low ) 



At, (6.9b) 



where 7^ low and 3>^ low are defined as in (|6.3p and l/°)u n is the term in parentheses in 
(|6.9bp ; the last lines remain the same as in the respective methods (|4.2p and (|5.17p . The 
modified CGMs for the single- and two-component equations are (|4.14|) and (|5.20|) , 
The three cases of increasing numerical stiffness are chosen as, 



For (fL5l) : mildly stiff V = 4, P = 2.1 (jji = 5.08); 

stiffer V = 4, P = 1.94 {n = 5.01); 

stiffest Vb = 6, P = 0.92 = 7.93); 



(6.10a) 



For (EI 



mildly stiff 
stiffer 
stiffest 



V = 4, = 1.50 (/i^ 1 ) = 5.10), P( 2 ) = 1.00 (/^( 2 ) = 5.92); 

F = 4, P( x ) = 0.50 OuW = 4.98), p( 2 ) = 1.50 (/i( 2 ) = 6.62); 

y = 6, P( 1 ) = 0.49 (/x^) = 7.93), P( 2 ) = 0.60 (/^ = 8.55) , 

(6.10b) 

where the values of //, as computed within the methods, are listed for reference only. Note 
that for the single-component equation, we had to choose the parameters for which the 
solution would satisfy the stability condition (|4.5p (see Fig. [T]) of the iterative methods 
(|4.2p . (|6.9ap . and ()4.14|) . For the two-component equation, the parameters were chosen by 
trial and error. 

The preconditioning operator was taken to be of the form (|2.18p or, in the two-component 



case, as 



N 



(c- V 2 ) diag(l, 1) 



(6.11) 



with c = 1 (which did not necessarily result in the optimal convergence rate) . The "fictitious 
time" step At for the ME-based methods was chosen as At = Ar op t — 0.1, as described 
above, where Ar op t was the optimal step size for the non-accelerated ITEM. The computa- 
tional domain and the initial conditions were taken as for the methods that find solutions 
with prescribed values of the propagation constant. Note that here, the modified CGMs 
(|4.14p and (|5.20p are guaranteed to converge to the solution for any initial conditions that 



resemble a two-dimensional pulse. Nonetheless we used ()6.8p with e x = 0.1 and e y 



-0.2 
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Iterative method 




non-accelerated (|4.2p 






Case in (|6.10a[) 


accel. by ME (j6.9a|) 


accel. by CGM (|4.14l) 


mildly stiff, At = 0.9 


330 iterations, 140 s 


90 iterations, 50 s 


50 iteratons, 40 s 


stiffer, At = 1.0 


1670 iterations, 710 s 


160 iterations, 100 s 


120 iteratons, 90 s 


stiffest, At = 0.6 


4690 iterations, 1980 s 


550 iterations, 330 s 


210 iteratons, 150 s 


Table 3: Comparative performance of the ITEM accelerated bv ME (16.9a[) and the modified 
CGM (|4.14P for the single-component Eq. (|1.5p. 




Iterative method 


Case in (l6.10bl) 


non-accelerated (|5.17p 


accel. by ME (l6.9b|) 


accel. by CGM (15.201) 


mildly stiff, At = 0.6 


300 iterations, 260 s 


120 iterations, 150 s 


60 iteratons, 80 s 


stiffer, At = 0.6 


850 iterations, 740 s 


220 iterations, 280 s 


120 iteratons, 130 s 


stiffest, At = 0.5 


1610 iterations, 1400 s 


380 iterations, 480 s 


130 iteratons, 160 s 



Table 4: Comparative performance of the ITEM accelerated by ME (|6.9bp and the modified 
CGM (j5T20T) for the two-component Eq. (jBTTj) . 



just for uniformity of all our simulations. For the same reason, we also started the acceler- 
ation at the same threshold (16. 6p . even though both the ME-based methods (16. 9p and the 
modified CGMs (|4,14p and (|5.20p could be employed at the very first iteration. The results 
are summarized in Tables 3 and 4. 

Again, as expected, we see that both ME- and CGM-based accelerations considerably 
improve the performance of the iterative method, with the stiffer the case, the more the 
improvement. Furthermore, the CGM-based acceleration performs better than the ME- 
based one by a factor of 2-3 for the two-component equation, with the stiffer the case, the 
more imrovement being provided by the CGM over the ME. Interestingly, however, for the 
single-component equation, the modified CGM provides an improvement over the ME only 
for the stiffest case; for somewhat less stiff cases, the improvement (in terms of computing 
time) is only marginal. 

7 Conclusions 

In this work, we proposed modifications of the well-known Conjugate Gradient method 
(CGM) that are capable of finding fundamental solitary waves of single- and multi-component 
Hamiltonian nonlinear equations. Our modified CGMs converge much faster than previously 
considered ierative methods of Richardson's type, like the (generalized) Petviashvili and 
imaginary-time evolution methods. Moreover, the slower the convergence of the Richardson- 
type method, the more speedup the modified CGM provides. 
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The classic CGM for a linear system of equations is known to converge when the matrix 
of this system is sign definite. For most solitary waves, the linearized operator about them, 
which is a counterpart of the matrix in the previous sentence, is not sign definite. While, 
in theory, this does not automatically imply that a straightforward generalization of the 
CGM to stationary nonlinear wave equations should fail (e.g., diverge), we verified that 
for the equations considered in Section 6, it does indeed fail. Then, the thrust of this 
work was to modify the CGM in such a way that it would be guaranteed to converge even 
when the linearized operator has eigenvalues of either sign. More precisely, our goal was to 
develop such modified versions of the CGM in the subcase when the number of the positive 
eigenvalues of that operator does not exceed the number of the components of the solitary 
wave. According to our "definition" at the end of Introduction, this situation would hold 
for fundamental solitary waves. 

The modified CGMs that we proposed in Sections 3 to 5 do not, strictly speaking, meet 
that goal. However, we show below that they come as closely as theoretically possible to 
doing so. As far as their practical applications, we demonstrated that our modified CGMs 
converge in the same parameter regions as earlier proposed, slower iterative methods. This 
convergence can be achieved by a simple shift of the initial condition, which has a theoretical 
explanation given in Section 6. 

When finding solitary waves with prescribed values of the propagation constants, we 
consider an equivalent nonlinear equation whose linearized operator is modified in such a 
way that its only positive eigenvalue is essentially turned into a negative one; see Eqs. (|3. 18|) 
and (|5.8p . This, however, makes the modified linearized operator "slightly" non-self-adjoint, 
because the generalized eigenfunction(s) of the original linearized operator, as in (|3.12p and 
(I5.4h . is (are) available only approximately. This is a fundamental fact about all nonlinear 
equations except for their very narrow subclass, equations with power-law nonlinearity 
(jl.ip (see also their two-component generalization in [7]), and hence it, in general, cannot 
be improved. Thus, this fact causes "slight" non-self- adjointness of the modified linearized 
operator. This, in turn, leads to sign indefiniteness of the quadratic form in (|3,17b .e) and 
thereby prevents one from obtaining a rigoroius guarantee that the modified CGM would 
always converge. However, we explained in Section 6 and Appendix 1 that the method can 
diverge only when L has small eigenvalues, and for those cases pointed out that a simple 
shift of the initial condition will suffice to make the modified CGM to actually converge. 

For finding solitary waves with prescribed values of the power (|1.7|) or, more generally, 
quadratic conserved quantities (|5.9p . the situation is different. There, the modified CGM 
is caused to converge not by modifying a linearized operator but by making the search 
directions satisfy a certain orthogonality relation; see (|4.3p and (|5.19p . The linearized 
operator employed by the method can be shown [13] to be self-adjoint in the space of 
functions satisfying those orthogonality relations, but it is negative definite only under 
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conditions stated before Eqs. (|4.5p and (|5.18p . This restriction is intrinsic to a method 
seeking a solution with a prescribed power, and hence cannot be relaxed for any nonlinear 
wave equations. 

Thus, we have justified why the modified CGMs proposed in this work come as closely 
as theoretically possible to guaranteeing that they would converge to fundamental solitary 
waves. In practice, however, these methods can always be forced to converge, as have been 
mentioned above and demonstrated in Section 6. 

The modified CGMs are faster not only than Richardson-type methods, but also than 
those methods accelerated by the slowest-decaying mode elimination (ME) technique [HE]. 
Namely, in comparison to the ME-accelerated methods, the modified CGMs provide an 
improvement of about a factor of two to three in terms of computing time; see Tables 
1-4 in Section 6 for details. Importantly, the CGMs do so in the most numerically stiff 
cases, when the respective non-accelerated methods would converge extremely slowly and 
hence the acceleration would be most desirable. In such cases, it would pay off to use the 
modified CGMs instead of ME, even though the latter is a little easier to program into a 
code. On the other hand, in non-stiff cases, i.e. when the non-accelerated methods would 
converge in just a few tens of iterations, both the modified CGMs and ME would provide 
only modest improvement in computing time; compare (|2.15p and (|2.2ip . In such cases, it 
may be simpler to use ME or even the non-accelerated method. Let us point out one other, 
practical, advantage of the ME-based acceleration over the CGMs. Namely, the fact that 
we were able to construct modified CGMs is critically grounded in the existence of relations 
(|3. 18|) . (|5.8p or (|4.3p . (|5.19p . as we explained above. On the contrary, ME can be used to 
accelerate any Richardson- type iterative method, e.g., methods from [3]-[6] or the ITEM 
with amplitude normalization |13j . Of course, while these methods do converge in many 
cases, their convergence conditions are not known, and hence their use with or without 
ME-based acceleration would come without a guarantee that they would converge. 

Finally, let us briefly mention alternative iterative methods, other than Newton's method 
proper, which can be used to compute solitary waves, both fundamental and non-fundamental 
In a recent numerical study |29j . Yang showed that a certain combination of the CGM and 
Newton's method converged to both fundamental and non-fundamental solitary waves for 
all of the examples considered in that study. It is remarkable, and at the moment has no 
rigorous theoretical explanation, that the straightforward generalization of the CGM out- 
lined at the end of Section 3.1 would diverge for most of the same examples for which Yang's 
method converges. Thus, Yang's method can be used in practice; however, it should be kept 
in mind that it can be possible to encounter situations where it would diverge. Alterna- 
tively, if one seeks a non-fundamental solitary wave, one can "square" the equation (see (|2,8p 
and (|2.9p ) and apply the original CGM to it. For linear systems, this trick has long been 
known (see, e.g., [19J, p. 304), and some researchers used it for finding non- fundamental 
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solitary waves [30]. Let us note, however, that "squaring" an equation leads to squaring the 
condition number of the corresponding linearized operator (since cond(A 2 ) = (cond(A)) 2 ), 
which, according to (|2.2ip . slows down the convergence; also, more arithmetic operatons 
per iteration are required for a "squared" method. 
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Appendix 1: Technical results for Sections 3 and 4 

Here we will first prove statements (i)-(iii) found in the first paragraph of Section 3.3 and 
also explain the choice (j3.15b|) for the constant T. Then we will prove that the fact that 
L may have a zero eigenvalue due to a translational eigenmode will not cause the CGM to 
break down. Finally, we will supply missing steps in the derivations of (|4.12p and (|4,13p . 
From (I3.15a|) . 

K (0,mo d ) v = ^ K (0) v = XVj (v,KWv) 

(v,v) 

Substituting the second equation in (jAl.ip into the definition of x we find: 

(v 1 Xv) 

X = F-, r = ^r. (A1.2) 

(v,v) 

Since T ^ 1 (see (j3. 15b|) ) . Eq. (|A1.2p implies that x = 0, which along with the second 
equation in (jAl.ip shows that (K^v = 0) 44> (K ( - > mod ^v = 0). This proves statement (i). 

Statement (ii), i.e. that i^( mod ) is approximately self-adjoint, follows by considering the 
following inner product for arbitrary functions / and g: 



\ I (v,v) 

» (f,Kg)-T 

» (g,Kf)-r 



(v,Kg) 
v, v) 

M^g) (vj) 

(v,v) 
(v,g) (v,Kf) 



(v,v) 

g,K^f) . (A1.3) 



The two approximate equalities above are due to the approximate relation (|3.12p . and we 
have also used the fact that K is self-adjoint. 

Statement (hi), i.e. that all eigenvalues of i^( mod ) that are not too close to zero are 
negative, follows by analogy with the statement found after Eq. (|2.13bp . First, due to 
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(|3.12p . one eigenfunction, V > of K^ mod ^ is close to v, and for it, 

K (mod)^(l) ^ ^(mod)^ _ A (l)( 1 _ r ^ _ A (l)( 1 _ r ^(l) _ (A1.4) 

When r is chosen according to (|3.15bp . the eigenvalue corresponding to Tp 1 ) is approxi- 
mately — 1, which is near the accumulation point of the spectrum of K |13j . It is convenient 
to "place" this eigenvalue inside (as opposed to at the edge of) the spectrum of iC( mod ) 
because then it does not affect the condition number of the operator and hence the conver- 
gence factor of the iterative method; see (]2.16p . (|2.15p . (12.21[) . Thus, we have shown that 
the "main culprit" of sign indefiniteness of K — the eigenvalue — has been successfully 
dealt with. To finish the proof of statement (iii), let us show that K ( mod ) can, in principle, 
acquire small positive eigenvalue(s). Let tft( k \ k > 2 be eigenfunctions of K with negative 
eigenvalues. Since v is only an approximate eigenfunction of the self-adjoint operator K, 
then (v,KiftW) = (Kv,ip( k ^) does not exactly equal to zero, which causes the eigenfunc- 
tions, and hence eigenvalues, of K^ mod ^ to differ slightly from those of K. Hence small 
negative eigenvalues of K could, in theory, become small positive eigenvalues of K^ mod \ 
This, however, does not appear to be the case in practice, since if it had been, then the 
generalized Petviashvili method (|3.7p would diverge (see the text before Eq. (|3.6p ). whereas 
our numerical experiments conducted in [7J and in this paper do not show such a divergence. 

Even if #( mod ) has only negative eigenvalues, the quadratic form (d n , K^ mod ^d n ) may 
still not be negative definite because i^( mod ) is not self-adjoint. Since, however, x( mod ) 
is only "slightly" non-self-adjoint, the sign indefiniteness of the above quadratic form can 
occur only when some eigenvalues of K^ mod ^ are close to zero, as illustrated by the following 
2x2 example: 

^'("i -i)(TH>°- (al5> 

Note that, for e <C 1, the vector in (|A1.5P is aligned primarily with the eigenvector, (5e 
corresponding to the smaller eigenvalue of the matrix. The practical implication of the sign 
indefiniteness of the quadratic form (d n , K^ mod ^d n ) is that it, and hence the denominators 
in (|3.17b .e). can become arbitrarily close to zero. This may cause divergence of the modified 
CGM ([XT?]) . 

To conclude the consideration of statement (iii), note that the eigenvalues of K = 
are the same as those of N l L. According to the Sylvester law of inertia 
|25j . eigenvalues of N~ l L and L have the same signs. Moreover, for reasonable choices of 
N (i.e., for c = O(l) in (I2.18|) ). the eigenvalues of both operators are of the same order 
of magnitude. Thus, small eigenvalues of K (and hence of K^ mod \ see above) imply small 
eigenvalues of L. 

Let us now show that if operator L has an eigenmode V'trans corresponding to transla- 
tional invariance of the solitary wave, this would not cause a breakdown of algorithm (|3.17p . 
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As before, we will perform the analysis in transformed variables (|3.1ip and f)3. 16j) and in 
the linear approximation with respect to v n , f n , and d n . 

The accordingly transformed translational eigenmode satisfies 

^ (m ° d) ^trans = 0. (A1.6) 

Let us also write down the counterparts, respectively numbered, of the lines of algorithm 
(|3.1|) that we will refer to: 

f n+ i = K^ mo ^v n+1 = K^v n+1 . (A1.7d) 



d n+1 = r n+1 + [5 n d n . (A1.7f) 
{f n+1 ,d n+ i) 

= - _ {K i^>i M i n+i) - _ < AL7b > 
We will show by contradiction that if d n is not proportional to Varans then neither will be 
dn+1. Then, clearly, since all the eigenvalues of Jf( mod ) other than V'trans are negative, the 
denominator in (|A1.7bjl will never vanish and hence the CGM will not break down. First, 
the solvability condition (Fredholm alternative) for (|A1.7dp implies that 

r n+ i = V?±, where (V>±, V'trans) = 0, (A1.8) 

and hence (3 n ^ since K^ mod ^d n ^ 0. Second, from §3?2§ and (1AT8J) it follows that 

{dn,$±)=0. (A1.9) 

Third, let us assume that for some n, d n+1 = a^ trans . Note that by (J3T2]) and (|A1.7f| . a ^ 0. 
Then from ()AL7fp it follows that 

dn = ^ tran ; " ^ => (dnM = ^±M^0. (ALIO) 

Pn Pn 

This contradicts (|A1.9h . and hence d n+ \ can never become proportional to V'trans, which 
proves our assertion. Note that if v n+ \ becomes proportional to V'trans , then f n+ \ = 0, in 
which case the iterations converge to (v + const • V'trans), which is just a translated solitary 
wave v. 

Finally, let us show that (|4,12p and (|4,13p hold. Since, as we noted in Section 4.2, the 
second equation in (|4.9cj) does not change the linearized form of the first equation, we can 
use v n +i instead of v n +i in (j4.9d|) . Then 

f n +i = r n + a n ICd n + 0(vl). (ALU) 

Along with (|4.9bp . this implies (|4.12p . Next, the only step that requires some explanation 
in the derivation of (|4.13p is the omission of (d n , N~ 1 v n +i). By (|4.10p . this inner product 
is 

(d n , iV~ V+i) = (4, ^~ V) + 0{v 2 n ) = 0{v 2 n ) , (A1.12) 
so that its omission does not change the accuracy implied by the right-hand side of (|4. 13|) . 
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Appendix 2: An alternative modified CGM for Section 3 



Since, for the reasons explained at the end of Section 3, we do not employ this algorithm 
for numerical examples presented in this paper, below we will present it only in terms of 
the transformed variables defined in (|3,lip and (|3,16p , 

f = K^v , d = f - ^iv , (A2.1a) 



(vo,vo) 
(v n ,Kd n ){r 



a. = . A ''' ( °" "'' > ■ < A2 - lb > 

(d n , Kd n ) 

- (vn, K«»v n ) 

v n+ i =v n + a n d n - —7-t- -v n , (A2.1c) 

A (1) \V n , V n ) 

K (o) 

v n+ \ (A2.1d) 

(f n+1 , K d n) _ (^fn)(f n+ uv n+1 ) 

n = ..- ^ Vn+l) , (A2.1e) 

T - | o 7 (v n +l, r n +l + Pnd n ) CAom 

dn+1 = r n +l + Pnd n ; : V n +1 ■ (A2.1f) 

{Vn+l,V n +l) 

The definitions of do m (|A2.1aP and d n+ i in ()A2.1fp ensure that 

(d n ,v n ) = (A2.2) 
at every iteration. Equation (|A2.1cP ensures that v n = v + v n where 

(v n ,v) — > as n increases . (A2.3) 
Equations (|A2.1b -d) entail the counterpart of (|3.2|) : 

(f n+ i, d n ) = 0[vl), (A2.4) 
and Eqs. (|A2.1e .f) entail a counterpart of (|3.3p : 

{d n+l , Kd n ) = 0{vl) . (A2.5) 
The last three relations can be proved similarly to relations (|4.12p and (|4.13p . 

Appendix 3: Sample code of modified CGM for a two-component 
solitary wave with prescribed propagation constants 

We will first present the code which uses the modified CGM (|3.17p . (|5.8p to find a solitary 
wave of Eqs. (|6.1|) . then explain some of its steps, and, finally, discuss how this code could 
be simplified. This code can be downloaded from 

http : //www. cems.uvm.edu/~lakobati/posted_papers_and_codes/ code_in_Appendix3 .m. 
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% Spatial and spectral domains: 

xlength=12*pi ; ylength=12*pi ; % domain lengths along x and y 

Nx=2~8; Ny=Nx; % number of points along x and y 

dx=xlength/Nx; dy=dx; % mesh sizes along x and y 

x=[-xlength/2:dx:xlength/2-dx] ; y=x; % domains along x and y 

[X,Y] =meshgrid(x,y) ; X and Y arrays of size Ny-by-Nx 

kx=2*pi/xlength* [0:Nx/2-l -Nx/2:-l]; ky=kx; % spectral domains along x and y 

[KX,KY] =meshgrid(kx,ky) ; '/, KX and KY arrays of size Ny-by-Nx 

DEL(: , : ,1)=-(KX.~2+1*KY.~2) ; DEL( : , : ,2)=DEL( : , : , 1) ; % Fourier symbol of Laplacian 

% Coefficients in the equation: 

mu(l)=7.89; mu(2)=8.5; % mu values 

W(: , : ,1)=6*( (cos(X)).~2 + (cos(Y)).~2 ) - mu(l) ; 

W(:,:,2)=6*( (cos(X)).~2 + (cos(Y)).~2 ) - mu(2) ; % potential minus mu 
F(l)=l; F12=0.5; F(2)=4; '/, nonlinearity coefficients 

Dt=0.9; 7, Delta tau 

% Initial condition: 

u(:,:,l)= 0.8*exp(-l*(X.~2 + Y. "2) ) . * (1+0 . l*X-0 . 2*Y) ; ; 
u(:,:,2)= 1.5*exp(-l*(X.~2 + Y. ~2) ) . * (1+0 . l*X-0 . 2*Y) ; ; 

% Loop control variables : 

normDu=l; % initialize the error norm to start the loop 

normDu_accel=0 . 05 ; % Acceleration begins when the error reaches this threshold; 

% parameters c, b, gamma are computed until this threshold. 
accelerate=0; % this marker indicates that the acceleration has started 

counter=0; °/ counter of the number of iterations 

% Iterate until the error reaches the prescribed tolerance 

while normDu >= 10" (-10) 
counter=counter+l ; 

if normDu >= normDu_accel & accelerate == °/ compute N, E2, gammas 

7 STEP 1: Compute parameters of N: b(2) , c(l) , c(2) 
DELu=real( ifft2(DEL.*fft2(u)) ); usq=u.~2; 
for k=l:2 

L0u(: , : ,k)=DELu(: , : ,k)+(W(: , : ,k)+F(k)*usq( : , : ,k)+F12*usq( : , : ,3-k)) .*u(: , : ,k) ; 
u_L0u(k)=trapz(trapz( u( : , : ,k) . *L0u( : , : ,k) )); % <u, L0u> componentwise 
end 

Lllul=L0u(: , : , l)+2*F(l)*usq( : , : ,1) .*u(: , : ,1) ; L12u2=2*F12*u( : , : ,1) .*usq(: , : ,2) ; 
L22u2=L0u(: , : , 2)+2*F(2) *usq( : , : ,2) .*u(: , : ,2); L21ul=2*F12*u( : , : ,2) .*usq(: , : ,1) ; 
SumjLkjEjK: , : , l)=Lllul+L12u2-L0u( : , : ,1) ; 
SumjLkjEjK: , : , 2)=L21ul+L22u2-L0u( : , : ,2) ; 
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for k=l:2 

u_SumjLkjEj 1 (k)=trapz(trapz( u( : , : ,k) . *SumjLkjEj 1 ( : , : ,k) )); 
DELu_Sum j LkjEj 1 (k) =trapz (tr apz ( DELu ( : , : , k) . *Sum j LkjEj 1 ( : , : ,k) ) ) ; 
u_u (k) =trapz (trapz ( u ( : , : , k) . *u ( : , : , k) )); 
u_DELu (k) =trapz (trapz ( u ( : , : , k) . *DELu ( : , : , k) ) ) ; 
DELu_DELu (k) =trapz (trapz ( DELu ( : , : , k) . *DELu ( : , : , k) ) ) ; 
end 

b(l)=l; 
for k=l:2 

kappa (k) = ( u_DELu (k) *DELu_Sum j Lk j E j 1 (k) -DELu_DELu (k) *u_Sum j Lk j E j 1 (k) )/. . . 
( u_u(k)*DELu_SumjLkjEjl(k)-u_DELu(k)*u_SumjLkjEjl(k) ); 

end 

b(2)=b(l)*(kappa(l)*u_u(l)-u_DELu(l)) * u_SumjLkjEj 1 (2) / ... 

( (kappa(2)*u_u(2)-u_DELu(2)) * u_SumjLkjEj 1 (1) ); 
c=b . *kappa; 

u_Nu=c.*u_u-b.*u_DELu; % <u, Nu> componentwise 

7, STEP 2: Compute eigenvector E2 = [rho2(l)*ul rho2 (2) *u2] "T and fft(N) ~ 
rho2(l)=-u_Nu(2)/u_Nu(l) ; rho2(2)=l; 

rhol(l)=l; rhol(2)=l; % these coefficients of El are for uniform notations 
for k=l:2 

El(: , : ,k)=rhol(k)* u (: , : ,k) ; E2(: , : ,k)=rho2(k)*u( : , : ,k) ; 
f f tN( : , : ,k)=c (k)-b(k) *DEL( : , : ,k) ; 7o Fourier symbol of N componentwise 
end 

7. STEP 3: Compute gammas 
E_NE(l)=sum( (rhol."2) .*u_Nu ); 
E_NE(2)=sum( (rho2 . "2) . * u _Nu ); 

SumjLkjEj2(: , : , l)=rho2(l) *Lllul+rho2(2) *L12u2-L0u( : , : ,1) ; 

SumjLkjEj2(: , : ,2)=rho2(l)*L21ul+rho2(2)*L22u2-L0u(: , : ,2) ; 

E_LE ( 1 ) =u_ Sum j Lk j E j 1 ( 1 ) +u_Sum j Lk j E j 1 ( 2 ) ; 

E_LE(2)=sum( trapz(trapz( E2 . *SumjLkjEj2 )) ); 

lambda=E_LE . /E_NE ; gamma=l+l . / (lambda*Dt) ; 

7 STEP 4: Update the solution: 

E_L0u ( 1 ) =rho 1 ( 1 ) *u_L0u ( 1) +rho 1 (2 ) *u_L0u (2 ) ; 

E_L0u (2) =rho2 ( 1 ) * u _L0u ( 1 ) +rho2 (2) *u_L0u (2) ; 

u = u + Dt*( real( if f t2 (f f t2 (LOu) . /f f tN) ) - ... 

gamma(l)*E_LOu(l)/E_NE(l)*El - gamma(2) *E_L0u(2) /E_NE(2) *E2 ); 

else 7o use last computed N, El ,E2, lambda, gamma and start CGM 

if accelerate == 
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accelerate = 1; 
for k=l:2 

L0u(: , : ,k)=real( if f t2(DEL( : , : ,k) . *f f t2(u( : , : ,k) ) ) ) + ... 

(W(: , : ,k)+F(k)*u(:, : ,k) . ~2+F12*u( : , : ,3-k) . ~2) .*u( : , : ,k) ; 
u_L0u (k) =trapz (trapz ( u( : , : , k) . *LOu ( : , : , k) ) ) ; 
GAMMA (k) =l+l/lambda(k) ; 
end 

E_L0u ( 1 ) =rho 1(1) *u_LOu ( 1 ) +rho 1(2) *u_LOu (2 ) ; 
E_LOu (2) =rho2 ( 1 ) *u_LOu ( 1 ) +rho2 (2) *u_LOu (2) ; 
NEl=real( if ft2(f ft2(El) . *f ftN) ); 
NE2=real( ifft2(fft2(E2) .*fftN) ); 

LL0u=L0u-GAMMA ( 1 ) *E_L0u ( 1 ) /E_NE ( 1 ) *NE1-GAMMA(2) *E_L0u (2) /E_NE (2) *NE2 ; 
r=real( ifft2(fft2(LL0u) . /fftN) ); d=r; 
end 

for k=l:2 

Ld(: , : ,k)=real( if f t2(DEL( : , : ,k) . *f f t2(d( : , : ,k) ) ) ) + ... 

(W(: , : ,k)+3*F(k)*u(: , : ,k) .~2+F12*u(: , : ,3-k) ."2) .*d(: , : ,k)+. . . 
2*F12*u(: , : ,k) .*u(: , : ,3-k) .*d(: , : ,3-k) ; 

end 

El_Ld=sum( trapz (trapz (El. *Ld)) ); E2_Ld=sum( trapz (trapz (E2 . *Ld) ) ); 
LLd=Ld-GAMMA ( 1 ) *El_Ld/E_NE ( 1 ) *NE1 -GAMMA (2 ) *E2_Ld/E_NE (2 ) *NE2 I 
Nr_d=sum( trapz (trapz (LLOu. *d) ) ); LLd_d=sum( trapz (trapz (d. *LLd) ) ); 
alpha = -Nr_d/LLd_d; 
u = u+alpha*d; 
for k=l:2 

L0u(: , : ,k)=real( if f t2(DEL( : , : ,k) . *f f t2 (u( : , : ,k) ) ) ) + ... 

(W( : , : ,k)+F(k)*u( : , : ,k) . "2+F12*u( : , : ,3-k) . ~2) . *u( : , : ,k) ; 
u_L0u (k) =trapz (trapz ( u ( : , : , k) . *LOu ( : , : , k) ) ) ; 
end 

E_L0u ( 1 ) =rho 1 ( 1 ) *u_L0u ( 1) +rho 1 (2 ) *u_LOu (2 ) ; 
E_L0u (2) =rho2 ( 1 ) *u_L0u ( 1 ) +rho2 (2) *u_L0u (2) ; 

LL0u=L0u-GAMMA ( 1 ) *E_L0u ( 1 ) /E_NE ( 1 ) *NE1-GAMMA (2) *E_L0u(2)/E_NE(2)*NE2 ; 
r = real( if ft2 (fft2 (LLOu) . /fftN) ); 
beta = max(-sum( trapz (trapz (r.*LLd)) )/LLd_d, 0); 
d = r + beta*d; 
end 

normDu= norm(L0u( : , : , 1) ) /norm(u( : , : , 1) )+norm(L0u( : , : , 2) ) /norm(u( : , : , 2) ) ; 
normDu_re corded (counter) =normDu ; 
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end 

figure(l); mesh(x,y ,u( : , : , 1) ) ; figure(2); mesh(x,y,u( : , : ,2) ) ; 
figure (3); plot ( [1 : counter] ,loglO(normDu_recorded)) 

Step 1 inside the while-loop implements Eqs. (4.12) and (4.13) (see also (4.21)) of [7J 
using the notations of that paper. For example, DELu, LOu, SumjLkjEj 1 in the code stand, re- 
spectively, for V 2 u n , Lou n , X]j=i Lkj&ji] note that all these quantities are two-dimensional 
vectors. Notations involving the underscore denote inner products; for example, u_L0u(k) 

(k) (k) 

denotes (uh, , (L(ju)n ), k = 1,2, where the superscript (k) means the feth component of 
the two-dimensional vector. 

Step 2 implements Eqs. (4.15)-(4.18) of [7J. Here rhol(k) denotes p^' k ' in the nota- 
tions of Eq. (|5.6p of this paper; note that, for the convenience of coding, the order of the 
superscripts is reversed compared to [TJ. 

Step 3 implements Eqs. (4.5) of [7J. Note that a& of [7] is denoted as in this 
paper. Also note that in the code, quantities like E_NE(1) denote (e^ 1 ), Ne' 1 '), i.e. here, 
the superscript refers to the particular eigenvector e^) of N _1 L. 

Step 4 implements Eq. (4.19) of [7j. 

At the first iteration of the CGM, in addition to implementing Eqs. (|3.17ap of this paper, 
one also computes and Ne^ fc \ k = 1,2, so that these quantities are not computed again 
at subsequent iterations of the CGM. Notations LLOu and LLd stand for L(°' mod )u n and 
L,( mod )d n . The remainder of the code implements Eqs. (|3.17p and (|5.8p of this paper. 

Note that the CGM iterations start when N, e( fc ), and are found rather imprecisely: 
accuracy (i.e., the value of normDu_accel) of 5% is used in the above code and in the 
examples reported in Section 6, and we verified that the code still worked when this accuracy 
was lowered to as much as 10%. Moreover, even if N etc. were computed up to a higher 
accuracy, the functions e^ k ' could still satisfy the eigenrelations (15. 4p only approximately. 
While this approximation is very close (99% in the least-squares sense) for e^, it is only 
about 70% for e( 2 ); see the second column in Table 1 of [7J. These considerations suggest 
that the code could still work if the N given by the expression above Eq. (|5.3p is replaced by 
a simpler expression (|6.1ip . The constant c there should be computed by a straightforward 
generalization of Eq. (3.11) of [7J, whereby the lines in Step 1 of the above code starting 
with the second for-loop through the end of that step are replaced by: 

for k=l:2 

u_u(k)=trapz(trapz( u( : , : ,k) . *u( : , : ,k) )); % <u, u> componentwise 

u_DELu(k)=trapz(trapz( u( : , : ,k) . *DELu( : , : ,k) )); % <u, DELu> componentwise 
end 

DELu_DELu=sum( trapz(trapz( DELu.*DELu )) ); % <DELu, DELu> as a scalar, etc 
u_SumjLkjEj l=sum( trapz(trapz( u. *SumjLkjEj 1 )) ); 
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DELu_SumjLkjEj l=sum( trapz(trapz( DELu. *SumjLkjEj 1 )) ); 

c = ( u_SumjLkjEjl*DELu_DELu - DELu_SumjLkjEj l*sum(u_DELu) )/... 

( u_SumjLkjEj l*sum(u_DELu) - DELu_SumjLkjEj l*sum(u_u) ); 
u_Nu=c*u_u-u_DELu; % <u, Nu> componentwise 

Note that the reason why the inner products u_u(k) and u_DELu(k) are computed for 
each component of u n rather than for the entire vector function, as DELu_DELu etc, is that 
they are needed to compute the individual components of u_Nu. The latter are needed to 
compute p( 2,1 ) in Step 2. 

We did not test the performance of such a simplified code because the goal of this paper 
is the proposal of modified CGMs for solitary waves and not optimization of those methods. 
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Figure 1: Schematics of the curves for fundamental solitary waves of Eq. (jl.5p with 

Vq = 4 and Vo = 6. The axes, markers, etc. are drawn not to scale. The three cases of 
increased numerical stiffness, listed in (I6.4ap . are labeled with filled circles. The shaded 
areas represent the first spectral bands of the linearized operator L for the two values of Vq. 
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Figure 3: Evolutions of the iteration error, denned by (|6.5p with 5 = 1, for the stiffest 
case (|6.4ap for Eq. (jl.5p . The curve for the non-accelerated generalized Petviashvili method 
(|3,7p is not shown in full so as not to obscure the details of the other two curves. 



44 



